import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import zscore
import seaborn as sns
sns.set_context('notebook', font_scale=1.5, rc={'lines.markeredgewidth': 2})
%load_ext autoreload
%autoreload 2
%matplotlib inline
import visual_behavior.visualization.utils as ut
import visual_behavior.data_access.loading as loading
import visual_behavior.visualization.ophys.summary_figures as sf
import visual_behavior.ophys.response_analysis.utilities as utilities
from visual_behavior.ophys.response_analysis.response_analysis import ResponseAnalysis
experiments_table = loading.get_filtered_ophys_experiment_table()
experiments_table = experiments_table[experiments_table.session_number.isin([2,5])==False]
cache_dir = loading.get_analysis_cache_dir()
df_name = 'stimulus_response_df'
conditions = ['cell_specimen_id', 'image_name', 'running', 'engagement_state']
project_codes = ['VisualBehaviorMultiscope']
experiments = experiments_table[experiments_table.project_code.isin(project_codes)]
stim_multi_session_df = loading.get_multi_session_df(cache_dir, df_name, conditions, experiments)
stim_multi_session_df.head()
cache_dir = loading.get_analysis_cache_dir()
df_name = 'trials_response_df'
conditions = ['cell_specimen_id', 'go' , 'change_image_name', 'running', 'engagement_state']
project_codes = ['VisualBehaviorMultiscope']
experiments = experiments_table[experiments_table.project_code.isin(project_codes)]
change_multi_session_df = loading.get_multi_session_df(cache_dir, df_name, conditions, experiments)
change_multi_session_df.head()
cache_dir = loading.get_analysis_cache_dir()
df_name = 'omission_response_df'
conditions = ['cell_specimen_id', 'running', 'engagement_state']
project_codes = ['VisualBehaviorMultiscope']
experiments = experiments_table[experiments_table.project_code.isin(project_codes)]
omission_multi_session_df = loading.get_multi_session_df(cache_dir, df_name, conditions, experiments)
omission_multi_session_df.head()
image_df = stim_multi_session_df.copy()
image_df = image_df[(image_df.running==True)&(image_df.engagement_state=='engaged')]
image_resp_cells = image_df[image_df.fraction_significant_p_value_gray_screen>0.2].cell_specimen_id.unique()
print(len(image_df), len(image_resp_cells))
omission_df = omission_multi_session_df.copy()
omission_df = omission_df[(omission_df.running==True)&(omission_df.engagement_state=='engaged')]
omission_resp_cells = omission_df[omission_df.fraction_significant_p_value_gray_screen>0.2].cell_specimen_id.unique()
print(len(omission_df), len(omission_resp_cells))
change_df = change_multi_session_df.copy()
change_df = change_df[(change_df.go==True)&(change_df.running==True)&(change_df.engagement_state=='engaged')]
change_resp_cells = change_df[change_df.fraction_significant_p_value_gray_screen>0.2].cell_specimen_id.unique()
print(len(change_df), len(change_resp_cells))
image_data = image_df[['cell_specimen_id', 'ophys_experiment_id', 'mean_trace', 'mean_response']]
image_data = image_data.rename(columns={'mean_trace':'image_trace', 'mean_response':'image_response'})
omission_data = omission_df[['cell_specimen_id', 'ophys_experiment_id', 'mean_trace', 'mean_response']]
omission_data = omission_data.rename(columns={'mean_trace':'omission_trace', 'mean_response':'omission_response'})
change_data = change_df[['cell_specimen_id', 'ophys_experiment_id', 'mean_trace', 'mean_response']]
change_data = change_data.rename(columns={'mean_trace':'change_trace', 'mean_response':'change_response'})
data = image_data.merge(omission_data, on=['cell_specimen_id', 'ophys_experiment_id'])
data = data.merge(change_data, on=['cell_specimen_id', 'ophys_experiment_id'])
data.shape
# add metadata to df
data = data.merge(experiments_table, on='ophys_experiment_id')
cells_to_include = np.unique(np.hstack((image_resp_cells, omission_resp_cells, change_resp_cells)))
print(len(stim_multi_session_df.cell_specimen_id.unique()), len(cells_to_include))
data = data[data.cell_specimen_id.isin(cells_to_include)]
data.shape
def get_cre_index(cre_line):
if cre_line == 'Slc17a7-IRES2-Cre':
cre_index = 0
elif cre_line == 'Sst-IRES-Cre':
cre_index = 1
elif cre_line == 'Vip-IRES-Cre':
cre_index = 2
return cre_index
data['cre_index'] = [get_cre_index(cre_line) for cre_line in data.cre_line.values]
# function to pull out just the traces matrix from dataframe
def get_concatenated_traces(data):
concatenated_traces = np.concatenate((np.vstack(data.image_trace.values), np.vstack(data.change_trace.values), np.vstack(data.omission_trace.values)), axis=1)
print(concatenated_traces.shape)
return concatenated_traces
cre_lines = np.sort(data.cre_line.unique())
len_images = image_df.mean_trace.values[0].shape[0]
len_changes= change_df.mean_trace.values[0].shape[0]
len_omissions = omission_df.mean_trace.values[0].shape[0]
fig, ax = plt.subplots(1,3, figsize=(15,5))
ax = ax.ravel()
for i,cre_line in enumerate(cre_lines):
order = np.argsort(data[data.cre_line==cre_line].image_response.values)
concatenated_traces = get_concatenated_traces(data[data.cre_line==cre_line])
concatenated_traces = concatenated_traces[order, :]
ax[i] = sns.heatmap(data=concatenated_traces, vmin=0, vmax=0.3, ax=ax[i], cbar=False)
ax[i].set_title(cre_line)
ax[i].axvline(x=len_images, ymin=0, ymax=1, color='white', linestyle='--')
ax[i].axvline(x=len_images+len_changes, ymin=0, ymax=1, color='white', linestyle='--')
ax[i].axvline(x=len_images+len_changes+len_omissions, ymin=0, ymax=1, color='white', linestyle='--')
fig.tight_layout()
fig, ax = plt.subplots(1,3, figsize=(15,5))
ax = ax.ravel()
for i,cre_line in enumerate(cre_lines):
order = np.argsort(data[data.cre_line==cre_line].image_response.values)
concatenated_traces = get_concatenated_traces(data[data.cre_line==cre_line])
concatenated_traces = concatenated_traces[order, :]
normed_traces = zscore(concatenated_traces, axis=1)
ax[i] = sns.heatmap(data=normed_traces, vmin=0, vmax=3, ax=ax[i], cbar=False)
ax[i].set_title(cre_line)
ax[i].axvline(x=len_images, ymin=0, ymax=1, color='white', linestyle='--')
ax[i].axvline(x=len_images+len_changes, ymin=0, ymax=1, color='white', linestyle='--')
fig.tight_layout()
import visual_behavior.clustering.multiscope_fn.def_funs as funs
concatenated_traces = get_concatenated_traces(data)
from sklearn.decomposition import PCA
varexpmax = .99 # 1 # .9
pc_all_cre = {}
pca_variance_all_cre = {}
print(f'Running PCA on matrix size: {np.shape(concatenated_traces)}')
x_train_pc, pca = funs.doPCA(concatenated_traces, varexpmax=varexpmax, doplot=1)
pca_variance = pca.explained_variance_ratio_
# x_train_pc.shape
pc_all_cre = x_train_pc
pca_variance_all_cre = pca_variance
fig, ax = plt.subplots(1,2,figsize=(10,5))
ax = ax.ravel()
c_value = data.cre_index.values
vmin = 0
vmax = 3
from matplotlib.colors import ListedColormap
cmap = ListedColormap(sns.color_palette('hsv', 3).as_hex())
for i,dim in enumerate(np.arange(1,3)):
x = pc_all_cre[:, dim]
y = pc_all_cre[:, dim+1]
cax = ax[i].scatter(x, y, s=2, marker='o', c=c_value, cmap=cmap, vmin=vmin, vmax=vmax)
ax[i].set_xlabel('PC '+str(dim))
ax[i].set_ylabel('PC '+str(dim+1))
ax[i].set_title('0=Slc, 1=Sst, 2=Vip')
# ax[i].set_xlim(-5,20)
# ax[i].set_ylim(-5,20)
cbar = plt.colorbar(cax)
cbar.ax.get_yaxis().labelpad = 15
cbar.ax.set_ylabel('cre line index', rotation=270)
fig.tight_layout()
fig, ax = plt.subplots(1,2,figsize=(10,5))
ax = ax.ravel()
c_value = data.omission_response.values
vmin = np.percentile(c_value, 5)
vmax = np.percentile(c_value, 95)
for i,dim in enumerate(np.arange(1,3)):
x = pc_all_cre[:, dim]
y = pc_all_cre[:, dim+1]
cax = ax[i].scatter(x, y, s=2, marker='o', c=c_value, cmap='magma', vmin=vmin, vmax=vmax)
ax[i].set_xlabel('PC '+str(dim))
ax[i].set_ylabel('PC '+str(dim+1))
ax[i].set_title('0=Slc, 1=Sst, 2=Vip')
cbar = plt.colorbar(cax)
cbar.ax.get_yaxis().labelpad = 15
cbar.ax.set_ylabel('mean omission response', rotation=270)
fig.tight_layout()
concatenated_traces = get_concatenated_traces(data)
concatenated_traces = zscore(concatenated_traces, axis=1)
from sklearn.decomposition import PCA
varexpmax = .99 # 1 # .9
pc_all_cre = {}
pca_variance_all_cre = {}
print(f'Running PCA on matrix size: {np.shape(concatenated_traces)}')
x_train_pc, pca = funs.doPCA(concatenated_traces, varexpmax=varexpmax, doplot=1)
pca_variance = pca.explained_variance_ratio_
# x_train_pc.shape
pc_all_cre = x_train_pc
pca_variance_all_cre = pca_variance
fig, ax = plt.subplots(1,2,figsize=(10,5))
ax = ax.ravel()
c_value = data.cre_index.values
vmin = 0
vmax = 3
from matplotlib.colors import ListedColormap
cmap = ListedColormap(sns.color_palette('hsv', 3).as_hex())
for i,dim in enumerate(np.arange(1,3)):
x = pc_all_cre[:, dim]
y = pc_all_cre[:, dim+1]
cax = ax[i].scatter(x, y, s=2, marker='o', c=c_value, cmap=cmap, vmin=vmin, vmax=vmax)
ax[i].set_xlabel('PC '+str(dim))
ax[i].set_ylabel('PC '+str(dim+1))
ax[i].set_title('0=Slc, 1=Sst, 2=Vip')
# ax[i].set_xlim(-5,20)
# ax[i].set_ylim(-5,20)
cbar = plt.colorbar(cax)
cbar.ax.get_yaxis().labelpad = 15
cbar.ax.set_ylabel('cre line index', rotation=270)
fig.tight_layout()
fig, ax = plt.subplots(1,2,figsize=(10,5))
ax = ax.ravel()
c_value = data.omission_response.values
vmin = np.percentile(c_value, 5)
vmax = np.percentile(c_value, 95)
for i,dim in enumerate(np.arange(1,3)):
x = pc_all_cre[:, dim]
y = pc_all_cre[:, dim+1]
cax = ax[i].scatter(x, y, s=2, marker='o', c=c_value, cmap='magma', vmin=vmin, vmax=vmax)
ax[i].set_xlabel('PC '+str(dim))
ax[i].set_ylabel('PC '+str(dim+1))
ax[i].set_title('0=Slc, 1=Sst, 2=Vip')
cbar = plt.colorbar(cax)
cbar.ax.get_yaxis().labelpad = 15
cbar.ax.set_ylabel('mean omission response', rotation=270)
fig.tight_layout()
import umap
import seaborn as sns
from mpl_toolkits.mplot3d import Axes3D
concatenated_traces = get_concatenated_traces(data)
ncomp = 3 # 3 # number of umap components
embedding_all_cre = {}
print(f'Running UMAP')
min_dist=0
sp = 1
neigh = 7
embedding = umap.UMAP(spread=sp, n_neighbors=neigh, n_components=ncomp,
min_dist=min_dist, metric='euclidean').fit_transform(concatenated_traces)
print(f'embedding size: {embedding.shape}')
def plot_UMAP_2D(embedding, c_value, cmap, vmin, vmax, label=None, title=None):
"""plots 3 components against each other in pairs - embedding must be 3 dimensional or greater"""
fig, ax = plt.subplots(1,2,figsize=(10,5))
ax = ax.ravel()
for i in range(2):
x = embedding[:, i]
y = embedding[:, i+1]
cax = ax[i].scatter(x, y, s=2, marker='o', c=c_value, cmap=cmap, vmin=vmin, vmax=vmax)
ax[i].set_xlabel('UMAP component '+str(i))
ax[i].set_ylabel('UMAP component '+str(i+1))
ax[i].set_title(title)
cbar = plt.colorbar(cax)
cbar.ax.get_yaxis().labelpad = 15
cbar.ax.set_ylabel(label, rotation=270)
fig.tight_layout()
data['cre_index'] = [get_cre_index(cre_line) for cre_line in data.cre_line.values]
c_value = data.cre_index.values
vmin = 0
vmax = 3
from matplotlib.colors import ListedColormap
cmap = ListedColormap(sns.color_palette('hsv', 3).as_hex())
title = '0=Slc, 1=Sst, 2=Vip'
label = 'cre_line index'
plot_UMAP_2D(embedding, c_value, cmap, vmin, vmax, label, title)
data.session_number = [int(session_type[6]) for session_type in data.session_type.values]
c_value = data.session_number.values
vmin = 0
vmax = 6
from matplotlib.colors import ListedColormap
cmap = ListedColormap(sns.color_palette('hsv', 6).as_hex())
title = 'colored by session'
label = 'session_number'
plot_UMAP_2D(embedding, c_value, cmap, vmin, vmax, label, title)
c_value = data.omission_response.values
vmin = np.percentile(c_value, 5)
vmax = np.percentile(c_value, 95)
cmap = 'magma'
label = 'mean_omission_response'
title = 'colored_by_omission'
plot_UMAP_2D(embedding, c_value, cmap, vmin, vmax, label, title)
c_value = data.change_response.values
vmin = np.percentile(c_value, 5)
vmax = np.percentile(c_value, 95)
cmap = 'magma'
title = 'colored by change response'
label = 'mean change response'
plot_UMAP_2D(embedding, c_value, cmap, vmin, vmax, label, title)
c_value = data.image_response.values
vmin = np.percentile(c_value, 5)
vmax = np.percentile(c_value, 95)
cmap = 'magma'
title = 'colored by image response'
label = 'mean image response'
plot_UMAP_2D(embedding, c_value, cmap, vmin, vmax, label, title)
import umap
import seaborn as sns
from mpl_toolkits.mplot3d import Axes3D
concatenated_traces = get_concatenated_traces(data)
concatenated_traces = zscore(concatenated_traces, axis=1)
ncomp = 3 # 3 # number of umap components
embedding_all_cre = {}
print(f'Running UMAP')
min_dist=0
sp = 5
neigh = 10
metric = 'euclidean'
embedding = umap.UMAP(spread=sp, n_neighbors=neigh, n_components=ncomp,
min_dist=min_dist, metric=metric).fit_transform(concatenated_traces)
print(f'embedding size: {embedding.shape}')
data['cre_index'] = [get_cre_index(cre_line) for cre_line in data.cre_line.values]
c_value = data.cre_index.values
vmin = 0
vmax = 3
from matplotlib.colors import ListedColormap
cmap = ListedColormap(sns.color_palette('hsv', 3).as_hex())
title = '0=Slc, 1=Sst, 2=Vip'
label = 'cre_line index'
plot_UMAP_2D(embedding, c_value, cmap, vmin, vmax, label, title)
data.session_number = [int(session_type[6]) for session_type in data.session_type.values]
c_value = data.session_number.values
vmin = 0
vmax = 6
from matplotlib.colors import ListedColormap
cmap = ListedColormap(sns.color_palette('hsv', 6).as_hex())
title = 'colored by session'
label = 'session_number'
plot_UMAP_2D(embedding, c_value, cmap, vmin, vmax, label, title)
c_value = data.omission_response.values
vmin = np.percentile(c_value, 5)
vmax = np.percentile(c_value, 95)
cmap = 'magma'
label = 'mean_omission_response'
title = 'colored_by_omission'
plot_UMAP_2D(embedding, c_value, cmap, vmin, vmax, label, title)
c_value = data.change_response.values
vmin = np.percentile(c_value, 5)
vmax = np.percentile(c_value, 95)
cmap = 'magma'
title = 'colored by change response'
label = 'mean change response'
plot_UMAP_2D(embedding, c_value, cmap, vmin, vmax, label, title)
c_value = data.image_response.values
vmin = np.percentile(c_value, 5)
vmax = np.percentile(c_value, 95)
cmap = 'magma'
title = 'colored by image response'
label = 'mean image response'
plot_UMAP_2D(embedding, c_value, cmap, vmin, vmax, label, title)